classdef Integration < handle

    properties
    
    end
    
    methods (Static)
        
        function [out, nodes, weights] = compute_percentiles(fct_handle, percentiles_to_compute, w_lb, w_ub, num_nodes, varargin)
            [function_mat,nodes,weights,w_transformed] = ...
                tools.Integration.prepare_integration(fct_handle, w_lb, w_ub, num_nodes, varargin{:});
           
            [nodes_sorted,idx]=sort(nodes);
            cdf = cumsum(function_mat(idx)'.*weights(idx));
            % compute percentiles in terms of nodes on [-1,1]
            percentiles_nodes = interp1q(cdf, nodes_sorted, ...
                                         percentiles_to_compute(:));
            % transform back to wages
            out = w_transformed(percentiles_nodes);
        end
        
       
        function [integral,nodes,weights] = integral_legendre(fct_handle, lb, ub, num_nodes, varargin)
        % integrates fct_handle over [lb, ub] using
        % gauss-legendre quadrature with num_nodes points
        % before integrating, the interval is transformed to [-1,1]
        % and w is transformed to log(w)
        %
        % Args:
        %    fct_handle (handle): function handle of wage 'density'
        %        over which to integrate. Needs to be a function of
        %        wage :math:`w` only.
        %    lb (double): lower-bound of wages over which to
        %        integrate, may be a vector
        %    ub (double): upper-bound of wages over which to
        %        integrate, may be a vector
        %    num_nodes (int): number of Gauss-Legendre nodes
        %    nodes (double, optional): pre-computed nodes, will be used if
        %         same length as num_nodes
        %     weights (double, optional): pre-computed weights, will be used if
        %         same length as num_nodes
        %    
        % Note:
        %     The size of the output is the same as
        %     :code:`size(lb)` and :code:`size(ub)`. If a
        %     column vector should be returned, the bounds have to
        %     be column vectors.
            
            p = inputParser;
            p.StructExpand = false;
            
            positive_num = @(w) isnumeric(w) && (min(w)>0);
            valid_num_nodes = @(x) isnumeric(x) && x==floor(x) && (x>0);
            
            addRequired(p, 'fct_handle', @(x) isa(x, 'function_handle'))
            addRequired(p, 'lb' , positive_num)
            addRequired(p, 'ub' , positive_num)
            addRequired(p, 'num_nodes', valid_num_nodes)
            addOptional(p, 'nodes', [], @isnumeric)
            addOptional(p, 'weights', [], @isnumeric)
            addParameter(p, 'plot_fct', false, @islogical)
            
            parse(p, fct_handle, lb, ub, num_nodes, varargin{:});
            
            nodes_pre_computed = p.Results.nodes;
            weights_pre_computed = p.Results.weights;
            
            [num_rows_lb, num_col_lb] = size(lb);
            [num_rows_ub, num_col_ub] = size(ub);

            assert(num_rows_lb == num_rows_ub && num_col_lb == num_col_ub, ['Tools:' ...
                                'Integration:integral_legendre:' ...
                                'lb_ub_different_size'],['ub and lb ' ...
                   'have different size'])
            
            assert(min(ub>lb)==1,'Tools:Integration:integral_legendre:lb_greater_ub','ub < lb');
            
            % computing nodes and weights is costly. Only do if
            % num_nodes differs from default, otherwise use
            % pre-computed nodes and weights
            if num_nodes == length(nodes_pre_computed)
                nodes = nodes_pre_computed;
                weights = weights_pre_computed;
            else
                warning('computing nodes and weights')
                [nodes, weights] = tools.Integration.lgwt(num_nodes,-1,1);
            end

            theta = @(phi,lb,ub) ...
              tools.ChangeOfVariables.change_variable_minone_one_log(phi,lb,ub);
        
            dtheta_dphi = @(phi,lb,ub) ...
            tools.ChangeOfVariables.change_variable_minone_one_log_deriv(phi,lb,ub);
            
            % bring bounds and nodes into matrix form
            lb_mat = repmat(lb(:),1,num_nodes);
            ub_mat = repmat(ub(:),1,num_nodes);
            nodes_mat = repmat(nodes(:)',length(lb),1);
            
            fun_at_nodes = fct_handle(theta(nodes_mat,lb_mat, ...
                                            ub_mat)).*dtheta_dphi(nodes_mat,lb_mat,ub_mat);
            
            tmp = fun_at_nodes * weights(:);            
            integral = reshape(tmp,size(lb));
        end

        function [function_mat,nodes,weights,w_transformed] = prepare_integration(fct_handle, w_lb, w_ub, num_nodes, varargin)
        % prepares integration as well as computation of
        % percentiles of fct_handle over [w_lb, w_ub] using
        % gauss-legendre quadrature with num_nodes points
        % before integrating, the interval is transformed to [-1,1]
        % and w is transformed to log(w)
        %
        % Args:
        %    fct_handle (handle): function handle of wage 'density'
        %        over which to integrate. Needs to be a function of
        %        wage :math:`w` only.
        %    w_lb (double): lower-bound of wages over which to
        %        integrate, may be a vector
        %    w_ub (double): upper-bound of wages over which to
        %        integrate, may be a vector
        %    num_nodes (int): number of Gauss-Legendre nodes
        %    nodes (double, optional): pre-computed nodes, will be used if
        %         same length as num_nodes
        %     weights (double, optional): pre-computed weights, will be used if
        %         same length as num_nodes
        %    
        % Keyword Arguments:
        %     plot_fct (logical): whether to plot the function
        %         value at the nodes, default is :code:`false`
        %
        % Note:
        %     The size of the output is the same as
        %     :code:`size(w_lb)` and :code:`size(w_ub)`. If a
        %     column vector should be returned, the bounds have to
        %     be column vectors.
            
        % ----- Input parsing
            p = inputParser;
            p.StructExpand = false;
            
            positive_num = @(w) isnumeric(w) && (min(w)>0);
            valid_num_nodes = @(x) isnumeric(x) && x==floor(x) && (x>0);
            
            addRequired(p, 'fct_handle', @(x) isa(x, 'function_handle'))
            addRequired(p, 'w_lb' , positive_num)
            addRequired(p, 'w_ub' , positive_num)
            addRequired(p, 'num_nodes', valid_num_nodes)
            addOptional(p, 'nodes', [], @isnumeric)
            addOptional(p, 'weights', [], @isnumeric)
            addParameter(p, 'plot_fct', false, @islogical)
            
            parse(p, fct_handle, w_lb, w_ub, num_nodes, varargin{:});
            
            nodes_pre_computed = p.Results.nodes;
            weights_pre_computed = p.Results.weights;
            
            plot_fct = p.Results.plot_fct;

            [num_rows_lb, num_col_lb] = size(w_lb);
            [num_rows_ub, num_col_ub] = size(w_ub);

            assert(num_rows_lb == num_rows_ub && num_col_lb == num_col_ub, ['Tools:' ...
                                'Integration:integral_legendre:' ...
                                'w_lb_w_ub_different_size'],['w_ub and w_lb ' ...
                   'have different size'])
            
            % -----
            % computing nodes and weights is costly. Only do if
            % num_nodes differs from default, otherwise use
            % pre-computed nodes and weights
            if num_nodes == length(nodes_pre_computed)
                nodes = nodes_pre_computed;
                weights = weights_pre_computed;
            else
                warning('computing nodes and weights')
                [nodes, weights] = tools.Integration.lgwt(num_nodes,-1,1);
            end
            
            num_nodes = length(nodes);
            num_w = max(num_rows_lb,num_col_lb);
            
            % loop over w_lb, w_ub
            jacobian_vector = zeros(num_nodes*num_w,1);
            nodes_transformed_vector = zeros(num_nodes*num_w,1);
            
            for i=1:num_w
                assert(w_ub(i) > w_lb(i), 'Tools:Integration:integral_legendre:w_lb_greater_w_ub','w_ub < w_lb')
                
                w_lb_i = w_lb(i);
                w_ub_i = w_ub(i);
                
                w_transformed = @(log_w) ...
                    tools.ChangeOfVariables ...
                    .change_variable_minone_one_log(log_w, w_lb_i, w_ub_i);
                
                dw_dlog_w = @(log_w) ...
                    tools.ChangeOfVariables.change_variable_minone_one_log_deriv(log_w,w_lb_i, w_ub_i);
                
                row_idx_start = (i-1)*num_nodes + 1;
                row_idx_end = i*num_nodes;
                
                nodes_transformed_vector(row_idx_start:row_idx_end,1) = w_transformed(nodes);
                jacobian_vector(row_idx_start:row_idx_end,1) = dw_dlog_w(nodes);
            end
            
            
            function_at_nodes = fct_handle(nodes_transformed_vector).*jacobian_vector;
            function_mat = reshape(function_at_nodes,[num_nodes,length(w_ub)])';
            
            if plot_fct
                % only works for scalar w_lb, w_ub
                scatter(nodes,function_at_nodes)
            end
        end        
        
        function mass = mass_bins(fct_handle, bins, num_nodes_bins,varargin)
        % compute mass of fct_handle over bins. Output is
        % vector of masses of length(bins)-1
        %
        % Args:
        %    fct_handle (handle): function handle of 'density'
        %        over which to integrate. Needs to be a function of
        %        x :math:`x` only.
        %     bins (double): vector of wage bin bounds
        %         where bins(1) to bins(2) is the first bin
        %         and bins(end-1) to bins(end) is the
        %         last bin
        %     num_nodes_bins (int): number of nodes to use for
        %         Gauss-Legendre integration within each wage bin
        %     nodes (double, optional): pre-computed nodes, will be used if
        %         same length as num_nodes
        %     weights (double, optional): pre-computed weights, will be used if
        %         same length as num_nodes
        %
        % Note:
        %     The size of the output is the same as
        %     :code:`size(bins)`.    
        
            
            
            % ----- Input parsing
            p = inputParser;
            p.StructExpand = false;
            
            positive_num = @(w) isnumeric(w) && (min(w)>0);
            validbins = @(bins) positive_num(bins) && ...
                issorted(bins) && isvector(bins) && ...
                length(bins)>=2;
            valid_num_nodes = @(x) isnumeric(x) && x==floor(x) && (x>0);
            
            addRequired(p, 'fct_handle', @(x) isa(x, 'function_handle'))
            addRequired(p, 'bins' , validbins)
            addRequired(p, 'num_nodes_bins' , valid_num_nodes)
            addOptional(p, 'nodes', [], @isnumeric)
            addOptional(p, 'weights', [], @isnumeric)
            
            parse(p, fct_handle, bins, num_nodes_bins, varargin{:});
            num_nodes_bins = p.Results.num_nodes_bins;
            nodes = p.Results.nodes;
            weights = p.Results.weights;
            
            w_lb = bins(1:end-1);
            w_ub = bins(2:end);
            
            mass = tools.Integration.integral_legendre(fct_handle, ...
                                                       w_lb, w_ub, ...
                                                       num_nodes_bins, varargin{:});
        end
        
        function [x,w]=lgwt(N,a,b)

        % lgwt.m
        %
        % This script is for computing definite integrals using Legendre-Gauss 
        % Quadrature. Computes the Legendre-Gauss nodes and weights  on an interval
        % [a,b] with truncation order N
        %
        % Suppose you have a continuous function f(x) which is defined on [a,b]
        % which you can evaluate at any x in [a,b]. Simply evaluate it at all of
        % the values contained in the x vector to obtain a vector f. Then compute
        % the definite integral using sum(f.*w);
        %
        % Written by Greg von Winckel - 02/25/2004
            N=N-1;
            N1=N+1; N2=N+2;

            xu=linspace(-1,1,N1)';

            % Initial guess
            y=cos((2*(0:N)'+1)*pi/(2*N+2))+(0.27/N1)*sin(pi*xu*N/N2);

            % Legendre-Gauss Vandermonde Matrix
            L=zeros(N1,N2);

            % Derivative of LGVM
            Lp=zeros(N1,N2);

            % Compute the zeros of the N+1 Legendre Polynomial
            % using the recursion relation and the Newton-Raphson method

            y0=2;

            % Iterate until new points are uniformly within epsilon of old points
            while max(abs(y-y0))>eps
                
                
                L(:,1)=1;
                Lp(:,1)=0;
                
                L(:,2)=y;
                Lp(:,2)=1;
                
                for k=2:N1
                    L(:,k+1)=( (2*k-1)*y.*L(:,k)-(k-1)*L(:,k-1) )/k;
                end
                
                Lp=(N2)*( L(:,N1)-y.*L(:,N2) )./(1-y.^2);   
                
                y0=y;
                y=y0-L(:,N2)./Lp;
                
            end

            % Linear map from[-1,1] to [a,b]
            x=(a*(1-y)+b*(1+y))/2;      

            % Compute the weights
            w=(b-a)./((1-y.^2).*Lp.^2)*(N2/N1)^2;        
        end
        
    end
end